Explore OBIS database and get species occurrences for Halibut (Hippoglossus hippoglossus).
halibut <- robis::occurrence(scientificname = "Hippoglossus hippoglossus",
geometry="POLYGON((-64.84626027288095 48.91641792598304,-64.01129933538095 50.49356505646415,-68.82331105413095 50.11464395818477,-71.04254933538095 46.80876488134784,-64.84626027288095 48.91641792598304))"
)## Retrieved 228 records of approximately 228 (100%)
halibut dataframe to a sf objectThe OBIS database provide a lot of context on the records.
## [1] "rightsHolder" "country"
## [3] "date_year" "institutionID"
## [5] "scientificNameID" "year"
## [7] "scientificName" "dynamicProperties"
## [9] "dropped" "aphiaID"
## [11] "language" "phylumid"
## [13] "familyid" "catalogNumber"
## [15] "basisOfRecord" "terrestrial"
## [17] "superclass" "modified"
## [19] "maximumDepthInMeters" "id"
## [21] "day" "order"
## [23] "superclassid" "dataset_id"
## [25] "collectionCode" "speciesid"
## [27] "date_end" "occurrenceID"
## [29] "license" "date_start"
## [31] "month" "genus"
## [33] "eventDate" "brackish"
## [35] "scientificNameAuthorship" "absence"
## [37] "subfamily" "genusid"
## [39] "originalScientificName" "marine"
## [41] "minimumDepthInMeters" "subphylumid"
## [43] "subfamilyid" "institutionCode"
## [45] "countryCode" "date_mid"
## [47] "eventTime" "class"
## [49] "orderid" "datasetName"
## [51] "kingdom" "specificEpithet"
## [53] "classid" "phylum"
## [55] "species" "subphylum"
## [57] "family" "category"
## [59] "kingdomid" "node_id"
## [61] "individualCount" "fieldNumber"
## [63] "type" "rights"
## [65] "recordedBy" "datasetID"
## [67] "locality" "taxonRank"
## [69] "depth" "geometry"
We perform a visual check of the occurrences.
We can now build a representation of the occurrences per year.
If you have a large dataset with a high density of points, it is a good idea to summarize the informations using a grid of abundance counts.
# get number of years
years <- unique(sf_halibut$year)
# Create a reference grid
ref_grid <- raster(sf_halibut, res=0.05)
# Open empty stack
stack_count <- stack()
for(i in 1:length(years)){
# Filter for one year (years[i])
sf_halibut_yr <- filter(sf_halibut, year == years[i])
# Prepare the geometry
sp_halibut_yr <- as(st_geometry(sf_halibut_yr), "Spatial")
# Count points per cell
rs_count <- rasterize(sp_halibut_yr, ref_grid, fun = 'count')
# Add the layer to the stack
stack_count <- addLayer(stack_count,rs_count)
}
names(stack_count) <- yearsrobis::occurrence(scientificname = "Hippoglossus hippoglossus")ggplot2.Hovmöller plot
data("NOAA_df_1990", package = "STRbook")
Tmax <- filter(NOAA_df_1990, # subset the data
proc == "Tmax" & # only max temperature
month %in% 5:9 & # May to September
year == 1993) # year of 1993
lim_lat <- range(Tmax$lat) # latitude range
lat_axis <- seq(lim_lat[1], # latitude axis
lim_lat[2],
length=25)
Tmax$t <- Tmax$julian - 728049 # create a new time variable
Tmax_grid <- Tmax
dists <- abs(outer(Tmax$lat, lat_axis, "-"))
Tmax_grid$lat <- lat_axis[apply(dists, 1, which.min)]
Tmax_lat_Hov <- group_by(Tmax_grid, lat, t) %>%
summarise(z = mean(z))
library(STRbook)
Hovmoller_lat <- ggplot(Tmax_lat_Hov) + # take data
geom_tile(aes(x = lat, y = t, fill = z)) + # plot
fill_scale(name = "degF") + # add color scale
scale_y_reverse() + # rev y scale
ylab("Day number (days)") + # add y label
xlab("Latitude (degrees)") + # add x label
theme_bw() # change theme
Hovmoller_latTrends plots by stations
data("NOAA_df_1990", package = "STRbook")
Tmax <- filter(NOAA_df_1990, # subset the data
proc == "Tmax" & # only max temperature
month %in% 5:9 & # May to September
year == 1993) # year of 1993
Tmax_av <- group_by(Tmax, date) %>%
summarise(meanTmax = mean(z))
gTmaxav <-
ggplot() +
geom_line(data = Tmax,aes(x = date, y = z, group = id),
colour = "blue", alpha = 0.04) +
geom_line(data = Tmax_av, aes(x = date, y = meanTmax)) +
xlab("Month") + ylab("Maximum temperature (degF)") +
theme_bw()
gTmaxavspacetime classTo perform spatio-temporal variograms with gstat, we need to convert our data into one of the following different spacetime classes of object. To build these objects, the best way is to use the stConstruct function in the spacetime package.
## station_id pm10 dates lng lat
## 1 DEBB051 8.688 2001-06-08 14.05746 52.54136
## 2 DEBB051 7.354 2001-09-10 14.05746 52.54136
## 3 DEBB051 15.729 2001-10-10 14.05746 52.54136
## 4 DEBB051 26.312 2001-05-03 14.05746 52.54136
## 5 DEBB051 13.958 2001-10-28 14.05746 52.54136
## 6 DEBB051 39.188 2001-10-16 14.05746 52.54136
spaceTime_air <- stConstruct(raw_air, c("lng","lat"),"dates", crs = CRS("+init=epsg:4326"))
class(spaceTime_air)## [1] "STIDF"
## attr(,"package")
## [1] "spacetime"
Class of spacetime
# subset
raw_air_062005 <- subset(raw_air, dates > as.Date("2005-06-01") & dates < as.Date("2005-06-30"))
raw_air_062005 <- stConstruct(raw_air_062005, c("lng","lat"),"dates", crs = CRS("+init=epsg:4326"))
# Coerce the STI class to STF
STFDF_air_062005 <- as(raw_air_062005, "STFDF")
summary(STFDF_air_062005)## Object of class STFDF
## with Dimensions (s, t, attr): (46, 28, 2)
## [[Spatial:]]
## Object of class SpatialPoints
## Coordinates:
## min max
## lng 6.28107 14.78617
## lat 47.80847 54.92497
## Is projected: FALSE
## proj4string :
## [+init=epsg:4326 +proj=longlat +datum=WGS84 +no_defs +ellps=WGS84
## +towgs84=0,0,0]
## Number of points: 46
## [[Temporal:]]
## Index timeIndex
## Min. :2005-06-02 Min. : 1.00
## 1st Qu.:2005-06-08 1st Qu.: 7.75
## Median :2005-06-15 Median :14.50
## Mean :2005-06-15 Mean :14.50
## 3rd Qu.:2005-06-22 3rd Qu.:21.25
## Max. :2005-06-29 Max. :28.00
## [[Data attributes:]]
## station_id pm10
## DESH001: 28 Min. : 2.708
## DENI063: 28 1st Qu.:10.739
## DEUB038: 28 Median :14.191
## DEBE056: 28 Mean :15.100
## DEBE032: 28 3rd Qu.:18.583
## (Other):1108 Max. :53.292
## NA's : 40 NA's :40